home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Gold Collection / Software Vault - The Gold Collection (American Databankers) (1993).ISO / cdr53 / acmalg01.zip / ACM600.FOR < prev    next >
Text File  |  1993-01-01  |  33KB  |  1,139 lines
  1. C     ALGORITHM 600, COLLECTED ALGORITHMS FROM ACM.
  2. C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.9, NO. 2,
  3. C     JUN., 1983, P. 258-259.
  4.       SUBROUTINE QUINAT(N, X, Y, B, C, D, E, F)
  5. C
  6.       INTEGER N
  7.       REAL X(N), Y(N), B(N), C(N), D(N), E(N), F(N)
  8. C
  9. C
  10. C
  11. C     QUINAT COMPUTES THE COEFFICIENTS OF A QUINTIC NATURAL QUINTIC SPLI
  12. C     S(X) WITH KNOTS X(I) INTERPOLATING THERE TO GIVEN FUNCTION VALUES:
  13. C               S(X(I)) = Y(I)  FOR I = 1,2, ..., N.
  14. C     IN EACH INTERVAL (X(I),X(I+1)) THE SPLINE FUNCTION S(XX) IS A
  15. C     POLYNOMIAL OF FIFTH DEGREE:
  16. C     S(XX) = ((((F(I)*P+E(I))*P+D(I))*P+C(I))*P+B(I))*P+Y(I)    (*)
  17. C           = ((((-F(I)*Q+E(I+1))*Q-D(I+1))*Q+C(I+1))*Q-B(I+1))*Q+Y(I+1)
  18. C     WHERE  P = XX - X(I)  AND  Q = X(I+1) - XX.
  19. C     (NOTE THE FIRST SUBSCRIPT IN THE SECOND EXPRESSION.)
  20. C     THE DIFFERENT POLYNOMIALS ARE PIECED TOGETHER SO THAT S(X) AND
  21. C     ITS DERIVATIVES UP TO S"" ARE CONTINUOUS.
  22. C
  23. C        INPUT:
  24. C
  25. C     N          NUMBER OF DATA POINTS, (AT LEAST THREE, I.E. N > 2)
  26. C     X(1:N)     THE STRICTLY INCREASING OR DECREASING SEQUENCE OF
  27. C                KNOTS.  THE SPACING MUST BE SUCH THAT THE FIFTH POWER
  28. C                OF X(I+1) - X(I) CAN BE FORMED WITHOUT OVERFLOW OR
  29. C                UNDERFLOW OF EXPONENTS.
  30. C     Y(1:N)     THE PRESCRIBED FUNCTION VALUES AT THE KNOTS.
  31. C
  32. C        OUTPUT:
  33. C
  34. C     B,C,D,E,F  THE COMPUTED SPLINE COEFFICIENTS AS IN (*).
  35. C         (1:N)  SPECIFICALLY
  36. C                B(I) = S'(X(I)), C(I) = S"(X(I))/2, D(I) = S"'(X(I))/6,
  37. C                E(I) = S""(X(I))/24,  F(I) = S""'(X(I))/120.
  38. C                F(N) IS NEITHER USED NOR ALTERED.  THE FIVE ARRAYS
  39. C                B,C,D,E,F MUST ALWAYS BE DISTINCT.
  40. C
  41. C        OPTION:
  42. C
  43. C     IT IS POSSIBLE TO SPECIFY VALUES FOR THE FIRST AND SECOND
  44. C     DERIVATIVES OF THE SPLINE FUNCTION AT ARBITRARILY MANY KNOTS.
  45. C     THIS IS DONE BY RELAXING THE REQUIREMENT THAT THE SEQUENCE OF
  46. C     KNOTS BE STRICTLY INCREASING OR DECREASING.  SPECIFICALLY:
  47. C
  48. C     IF X(J) = X(J+1) THEN S(X(J)) = Y(J) AND S'(X(J)) = Y(J+1),
  49. C     IF X(J) = X(J+1) = X(J+2) THEN IN ADDITION S"(X(J)) = Y(J+2).
  50. C
  51. C     NOTE THAT S""(X) IS DISCONTINUOUS AT A DOUBLE KNOT AND, IN
  52. C     ADDITION, S"'(X) IS DISCONTINUOUS AT A TRIPLE KNOT.  THE
  53. C     SUBROUTINE ASSIGNS Y(I) TO Y(I+1) IN THESE CASES AND ALSO TO
  54. C     Y(I+2) AT A TRIPLE KNOT.  THE REPRESENTATION (*) REMAINS
  55. C     VALID IN EACH OPEN INTERVAL (X(I),X(I+1)).  AT A DOUBLE KNOT,
  56. C     X(J) = X(J+1), THE OUTPUT COEFFICIENTS HAVE THE FOLLOWING VALUES:
  57. C       Y(J) = S(X(J))          = Y(J+1)
  58. C       B(J) = S'(X(J))         = B(J+1)
  59. C       C(J) = S"(X(J))/2       = C(J+1)
  60. C       D(J) = S"'(X(J))/6      = D(J+1)
  61. C       E(J) = S""(X(J)-0)/24     E(J+1) = S""(X(J)+0)/24
  62. C       F(J) = S""'(X(J)-0)/120   F(J+1) = S""'(X(J)+0)/120
  63. C     AT A TRIPLE KNOT, X(J) = X(J+1) = X(J+2), THE OUTPUT
  64. C     COEFFICIENTS HAVE THE FOLLOWING VALUES:
  65. C       Y(J) = S(X(J))         = Y(J+1)    = Y(J+2)
  66. C       B(J) = S'(X(J))        = B(J+1)    = B(J+2)
  67. C       C(J) = S"(X(J))/2      = C(J+1)    = C(J+2)
  68. C       D(J) = S"'((X(J)-0)/6    D(J+1) = 0  D(J+2) = S"'(X(J)+0)/6
  69. C       E(J) = S""(X(J)-0)/24    E(J+1) = 0  E(J+2) = S""(X(J)+0)/24
  70. C       F(J) = S""'(X(J)-0)/120  F(J+1) = 0  F(J+2) = S""'(X(J)+0)/120
  71. C
  72.       INTEGER I, M
  73.       REAL B1, P, PQ, PQQR, PR, P2, P3, Q, QR, Q2, Q3, R, R2, S, T, U, V
  74. C
  75.       IF (N.LE.2) GO TO 190
  76. C
  77. C     COEFFICIENTS OF A POSITIVE DEFINITE, PENTADIAGONAL MATRIX,
  78. C     STORED IN D,E,F FROM 2 TO N-2.
  79. C
  80.       M = N - 2
  81.       Q = X(2) - X(1)
  82.       R = X(3) - X(2)
  83.       Q2 = Q*Q
  84.       R2 = R*R
  85.       QR = Q + R
  86.       D(1) = 0.
  87.       E(1) = 0.
  88.       D(2) = 0.
  89.       IF (Q.NE.0.) D(2) = 6.*Q*Q2/(QR*QR)
  90. C
  91.       IF (M.LT.2) GO TO 40
  92.       DO 30 I=2,M
  93.         P = Q
  94.         Q = R
  95.         R = X(I+2) - X(I+1)
  96.         P2 = Q2
  97.         Q2 = R2
  98.         R2 = R*R
  99.         PQ = QR
  100.         QR = Q + R
  101.         IF (Q) 20, 10, 20
  102.    10   D(I+1) = 0.
  103.         E(I) = 0.
  104.         F(I-1) = 0.
  105.         GO TO 30
  106.    20   Q3 = Q2*Q
  107.         PR = P*R
  108.         PQQR = PQ*QR
  109.         D(I+1) = 6.*Q3/(QR*QR)
  110.         D(I) = D(I) + (Q+Q)*(15.*PR*PR+(P+R)*Q*(20.*PR+7.*Q2)+Q2*(8.*
  111.      *   (P2+R2)+21.*PR+Q2+Q2))/(PQQR*PQQR)
  112.         D(I-1) = D(I-1) + 6.*Q3/(PQ*PQ)
  113.         E(I) = Q2*(P*QR+3.*PQ*(QR+R+R))/(PQQR*QR)
  114.         E(I-1) = E(I-1) + Q2*(R*PQ+3.*QR*(PQ+P+P))/(PQQR*PQ)
  115.         F(I-1) = Q3/PQQR
  116.    30 CONTINUE
  117. C
  118.    40 IF (R.NE.0.) D(M) = D(M) + 6.*R*R2/(QR*QR)
  119. C
  120. C     FIRST AND SECOND ORDER DIVIDED DIFFERENCES OF THE GIVEN FUNCTION
  121. C     VALUES, STORED IN B FROM 2 TO N AND IN C FROM 3 TO N
  122. C     RESPECTIVELY. CARE IS TAKEN OF DOUBLE AND TRIPLE KNOTS.
  123. C
  124.       DO 60 I=2,N
  125.         IF (X(I).NE.X(I-1)) GO TO 50
  126.         B(I) = Y(I)
  127.         Y(I) = Y(I-1)
  128.         GO TO 60
  129.    50   B(I) = (Y(I)-Y(I-1))/(X(I)-X(I-1))
  130.    60 CONTINUE
  131.       DO 80 I=3,N
  132.         IF (X(I).NE.X(I-2)) GO TO 70
  133.         C(I) = B(I)*0.5
  134.         B(I) = B(I-1)
  135.         GO TO 80
  136.    70   C(I) = (B(I)-B(I-1))/(X(I)-X(I-2))
  137.    80 CONTINUE
  138. C
  139. C     SOLVE THE LINEAR SYSTEM WITH C(I+2) - C(I+1) AS RIGHT-HAND SIDE.
  140. C
  141.       IF (M.LT.2) GO TO 100
  142.       P = 0.
  143.       C(1) = 0.
  144.       E(M) = 0.
  145.       F(1) = 0.
  146.       F(M-1) = 0.
  147.       F(M) = 0.
  148.       C(2) = C(4) - C(3)
  149.       D(2) = 1./D(2)
  150. C
  151.       IF (M.LT.3) GO TO 100
  152.       DO 90 I=3,M
  153.         Q = D(I-1)*E(I-1)
  154.         D(I) = 1./(D(I)-P*F(I-2)-Q*E(I-1))
  155.         E(I) = E(I) - Q*F(I-1)
  156.         C(I) = C(I+2) - C(I+1) - P*C(I-2) - Q*C(I-1)
  157.         P = D(I-1)*F(I-1)
  158.    90 CONTINUE
  159. C
  160.   100 I = N - 1
  161.       C(N-1) = 0.
  162.       C(N) = 0.
  163.       IF (N.LT.4) GO TO 120
  164.       DO 110 M=4,N
  165. C        I = N-2, ..., 2
  166.         I = I - 1
  167.         C(I) = (C(I)-E(I)*C(I+1)-F(I)*C(I+2))*D(I)
  168.   110 CONTINUE
  169. C
  170. C     INTEGRATE THE THIRD DERIVATIVE OF S(X).
  171. C
  172.   120 M = N - 1
  173.       Q = X(2) - X(1)
  174.       R = X(3) - X(2)
  175.       B1 = B(2)
  176.       Q3 = Q*Q*Q
  177.       QR = Q + R
  178.       IF (QR) 140, 130, 140
  179.   130 V = 0.
  180.       T = 0.
  181.       GO TO 150
  182.   140 V = C(2)/QR
  183.       T = V
  184.   150 F(1) = 0.
  185.       IF (Q.NE.0.) F(1) = V/Q
  186.       DO 180 I=2,M
  187.         P = Q
  188.         Q = R
  189.         R = 0.
  190.         IF (I.NE.M) R = X(I+2) - X(I+1)
  191.         P3 = Q3
  192.         Q3 = Q*Q*Q
  193.         PQ = QR
  194.         QR = Q + R
  195.         S = T
  196.         T = 0.
  197.         IF (QR.NE.0.) T = (C(I+1)-C(I))/QR
  198.         U = V
  199.         V = T - S
  200.         IF (PQ) 170, 160, 170
  201.   160   C(I) = C(I-1)
  202.         D(I) = 0.
  203.         E(I) = 0.
  204.         F(I) = 0.
  205.         GO TO 180
  206.   170   F(I) = F(I-1)
  207.         IF (Q.NE.0.) F(I) = V/Q
  208.         E(I) = 5.*S
  209.         D(I) = 10.*(C(I)-Q*S)
  210.         C(I) = D(I)*(P-Q) + (B(I+1)-B(I)+(U-E(I))*P3-(V+E(I))*Q3)/PQ
  211.         B(I) = (P*(B(I+1)-V*Q3)+Q*(B(I)-U*P3))/PQ -
  212.      *   P*Q*(D(I)+E(I)*(Q-P))
  213.   180 CONTINUE
  214. C
  215. C     END POINTS X(1) AND X(N).
  216. C
  217.       P = X(2) - X(1)
  218.       S = F(1)*P*P*P
  219.       E(1) = 0.
  220.       D(1) = 0.
  221.       C(1) = C(2) - 10.*S
  222.       B(1) = B1 - (C(1)+S)*P
  223. C
  224.       Q = X(N) - X(N-1)
  225.       T = F(N-1)*Q*Q*Q
  226.       E(N) = 0.
  227.       D(N) = 0.
  228.       C(N) = C(N-1) + 10.*T
  229.       B(N) = B(N) + (C(N)-T)*Q
  230.   190 RETURN
  231.       END
  232. C
  233. C
  234.       SUBROUTINE QUINEQ(N, Y, B, C, D, E, F)
  235. C
  236.       INTEGER N
  237.       REAL Y(N), B(N), C(N), D(N), E(N), F(N)
  238. C
  239. C
  240. C
  241. C     QUINEQ COMPUTES THE COEFFICIENTS OF QUINTIC NATURAL QUINTIC SPLINE
  242. C     S(X) WITH EQUIDISTANT KNOTS X(I) INTERPOLATING THERE TO GIVEN
  243. C     FUNCTION VALUES:
  244. C               S(X(I)) = Y(I)  FOR I = 1,2, ..., N.
  245. C     IN EACH INTERVAL (X(I),X(I+1)) THE SPLINE FUNCTION S(XX) IS
  246. C     A POLYNOMIAL OF FIFTH DEGREE:
  247. C     S(XX)=((((F(I)*P+E(I))*P+D(I))*P+C(I))*P+B(I))*P+Y(I)    (*)
  248. C          =((((-F(I)*Q+E(I+1))*Q-D(I+1))*Q+C(I+1))*Q-B(I+1))*Q+Y(I+1)
  249. C     WHERE  P = (XX - X(I))/X(I+1) - X(I))
  250. C     AND    Q = (X(I+1) - XX)/(X(I+1) - X(I)).
  251. C     (NOTE THE FIRST SUBSCRIPT IN THE SECOND EXPRESSION.)
  252. C     THE DIFFERENT POLYNOMIALS ARE PIECED TOGETHER SO THAT S(X) AND
  253. C     ITS DERIVATIVES UP TO S"" ARE CONTINUOUS.
  254. C
  255. C        INPUT:
  256. C
  257. C     N          NUMBER OF DATA POINTS, (AT LEAST THREE, I.E. N > 2)
  258. C     Y(1:N)     THE PRESCRIBED FUNCTION VALUES AT THE KNOTS
  259. C
  260. C        OUTPUT:
  261. C
  262. C     B,C,D,E,F  THE COMPUTED SPLINE COEFFICIENTS AS IN (*).
  263. C         (1:N)  IF X(I+1) - X(I) = 1., THEN SPECIFICALLY:
  264. C                B(I) = S'X(I)), C(I) = S"(X(I))/2, D(I) = S"'(X(I))/6,
  265. C                E(I) = S""(X(I))/24,  F(I) = S""'(X(I)+0)/120.
  266. C                F(N) IS NEITHER USED NOR ALTERED.  THE ARRAYS
  267. C                Y,B,C,D MUST ALWAYS BE DISTINCT.  IF E AND F ARE
  268. C                NOT WANTED, THE CALL QUINEQ(N,Y,B,C,D,D,D) MAY
  269. C                BE USED TO SAVE STORAGE LOCATIONS.
  270. C
  271.       INTEGER I, M
  272.       REAL P, Q, R, S, T, U, V
  273. C
  274.       IF (N.LE.2) GO TO 50
  275. C
  276.       M = N - 3
  277.       P = 0.
  278.       Q = 0.
  279.       R = 0.
  280.       S = 0.
  281.       T = 0.
  282.       D(M+1) = 0.
  283.       D(M+2) = 0.
  284.       IF (M.LE.0) GO TO 30
  285.       DO 10 I=1,M
  286.         U = P*R
  287.         B(I) = 1./(66.-U*R-Q)
  288.         R = 26. - U
  289.         C(I) = R
  290.         D(I) = Y(I+3) - 3.*(Y(I+2)-Y(I+1)) - Y(I) - U*S - Q*T
  291.         Q = P
  292.         P = B(I)
  293.         T = S
  294.         S = D(I)
  295.    10 CONTINUE
  296. C
  297.       I = N - 2
  298.       DO 20 M=4,N
  299. C        I    = N-3, ..., 1
  300.         I = I - 1
  301.         D(I) = (D(I)-C(I)*D(I+1)-D(I+2))*B(I)
  302.    20 CONTINUE
  303. C
  304.    30 M = N - 1
  305.       Q = 0.
  306.       R = D(1)
  307.       T = R
  308.       V = R
  309.       DO 40 I=2,M
  310.         P = Q
  311.         Q = R
  312.         R = D(I)
  313.         S = T
  314.         T = P - Q - Q + R
  315.         F(I) = T
  316.         U = 5.*(-P+Q)
  317.         E(I) = U
  318.         D(I) = 10.*(P+Q)
  319.         C(I) = 0.5*(Y(I+1)+Y(I-1)+S-T) - Y(I) - U
  320.         B(I) = 0.5*(Y(I+1)-Y(I-1)-S-T) - D(I)
  321.    40 CONTINUE
  322. C
  323.       F(1) = V
  324.       E(1) = 0.
  325.       D(1) = 0.
  326.       C(1) = C(2) - 10.*V
  327.       B(1) = Y(2) - Y(1) - C(1) - V
  328.       E(N) = 0.
  329.       D(N) = 0.
  330.       C(N) = C(N-1) + 10.*T
  331.       B(N) = Y(N) - Y(N-1) + C(N) - T
  332.    50 RETURN
  333.       END
  334. C
  335. C
  336. C
  337. C
  338.       SUBROUTINE QUINDF(N, X, Y, B, C, D, E, F)
  339. C
  340.       INTEGER N
  341.       REAL X(N), Y(N), B(N), C(N), D(N), E(N), F(N)
  342. C
  343. C
  344. C
  345. C     QUINDF COMPUTES THE COEFFICIENTS OF A QUINTIC NATURAL QUINTIC
  346. C     SPLINE S(X) WITH KNOTS X(I) FOR WHICH THE FUNCTION VALUES Y(I) AND
  347. C     THE FIRST DERIVATIVES B(I) ARE SPECIFIED AT X(I), I = 1,2, ...,N.
  348. C     IN EACH INTERVAL (X(I),X(I+1)) THE SPLINE FUNCTION S(XX) IS
  349. C     A POLYNOMIAL OF FIFTH DEGREE:
  350. C     S(XX)=((((F(I)*P+E(I))*P+D(I))*P+C(I))*P+B(I))*P+Y(I)    (*)
  351. C     WHERE  P = XX - X(I).
  352. C
  353. C        INPUT:
  354. C
  355. C     N          NUMBER OF DATA POINTS, (AT LEAST TWO, I.E. N > 1)
  356. C     X(1:N)     THE STRICTLY INCREASING OR DECREASING SEQUENCE OF
  357. C                KNOTS.  THE SPACING MUST BE SUCH THAT THE FIFTH
  358. C                POWER OF X(I+1) - X(I) CAN BE FORMED WITHOUT
  359. C                OVERFLOW OR UNDERFLOW OF EXPONENTS
  360. C     Y(1:N)     THE PRESCRIBED FUNCTION VALUES AT THE KNOTS
  361. C     B(1:N)     THE PRESCRIBED DERIVATIVE VALUES AT THE KNOTS
  362. C
  363. C        OUTPUT:
  364. C
  365. C     C,D,E,F    THE COMPUTED SPLINE COEFFICIENTS AS IN (*).
  366. C         (1:N)  E(N) AND F(N) ARG NEITHER USED NOR ALTERED.
  367. C                THE ARRAYS C,D,E,F MUST ALWAYS BE DISTINCT.
  368. C
  369.       INTEGER I, M, N1
  370.       REAL CC, G, H, HH, H2, P, PP, Q, QQ, R, RR
  371. C
  372.       IF (N.LE.1) GO TO 40
  373.       N1 = N - 1
  374.       CC = 0.
  375.       HH = 0.
  376.       PP = 0.
  377.       QQ = 0.
  378.       RR = 0.
  379.       G = 0.
  380.       DO 10 I=1,N1
  381.         H = 1./(X(I+1)-X(I))
  382.         H2 = H*H
  383.         D(I) = 3.*(HH+H) - G*HH
  384.         P = (Y(I+1)-Y(I))*H2*H
  385.         Q = (B(I+1)+B(I))*H2
  386.         R = (B(I+1)-B(I))*H2
  387.         CC = 10.*(P-PP) - 5.*(Q-QQ) + R + RR + G*CC
  388.         C(I) = CC
  389.         G = H/D(I)
  390.         HH = H
  391.         PP = P
  392.         QQ = Q
  393.         RR = R
  394.    10 CONTINUE
  395. C
  396.       C(N) = (-10.*PP+5.*QQ+RR+G*CC)/(3.*HH-G*HH)
  397.       I = N
  398.       DO 20 M=1,N1
  399. C        I      = N-1, ..., 1
  400.         I = I - 1
  401.         D(I+1) = 1./(X(I+1)-X(I))
  402.         C(I) = (C(I)+C(I+1)*D(I+1))/D(I)
  403.    20 CONTINUE
  404. C
  405.       DO 30 I=1,N1
  406.         H = D(I+1)
  407.         P = (((Y(I+1)-Y(I))*H-B(I))*H-C(I))*H
  408.         Q = ((B(I+1)-B(I))*H-C(I)-C(I))*H
  409.         R = (C(I+1)-C(I))*H
  410.         G = Q - 3.*P
  411.         RR = R - 3.*(P+G)
  412.         QQ = -RR - RR + G
  413.         F(I) = RR*H*H
  414.         E(I) = QQ*H
  415.         D(I) = -RR - QQ + P
  416.    30 CONTINUE
  417. C
  418.       D(N) = 0.
  419.       E(N) = 0.
  420.       F(N) = 0.
  421.    40 RETURN
  422.       END
  423. C    DRIVER PROGRAM FOR TEST OF QUINAT, QUINEQ AND QUINDF
  424. C       FOLLOWS.
  425. C
  426.       INTEGER N,NM1,M,MM,MM1,I,K,J,JJ
  427.       REAL Z
  428.       REAL X(200),Y(200),B(200),BB(200),CC(200),DD(200),EE(200),
  429.      *                 FF(200),A(200,6),C(6),DIFF(5),COM(5)
  430. C
  431. C     N          NUMBER OF DATA POINTS.
  432. C     M          2*M-1 IS ORDER OF SPLINE.
  433. C                   M = 3 ALWAYS FOR QUINTIC SPLINE.
  434. C     NN,NM1,MM,
  435. C     MM1,I,K,
  436. C     J,JJ       TEMPORARY INTEGER VARIABLES.
  437. C     Z,P        TEMPORARY REAL VARIABLES.
  438. C     X(1:N)     THE SEQUENCE OF KNOTS.
  439. C     Y(1:N)     THE PRESCRIBED FUNCTION VALUES AT THE KNOTS.
  440. C     B(1:N)     THE PRESCRIBED DERIVATIVE VALUES AT THE KNOTS.
  441. C     BB,CC,DD,
  442. C     EE,FF(1:N) THE COMPUTED SPLINE COEFFICIENTS
  443. C     A(1:N,1:6) TWO DIMENSIONAL ARRAY WHOSE COLUMNS ARE
  444. C                   Y, BB, CC, DD, EE, FF.
  445. C     DIFF(1:5)  MAXIMUM VALUES OF DIFFERENCES OF VALUES AND
  446. C                   DERIVATIVES TO RIGHT AND LEFT OF KNOTS.
  447. C     COM(1:5)   MAXIMUM VALUES OF COEFFICIENTS.
  448. C
  449. C
  450. C     TEST OF QUINAT WITH NONEQUIDISTANT KNOTS AND
  451. C        EQUIDISTANT KNOTS FOLLOWS.
  452. C
  453.       NOUT=6
  454.       WRITE(NOUT,10)
  455.    10 FORMAT(50H1         TEST OF QUINAT WITH NONEQUIDISTANT KNOTS)
  456.       N = 5
  457.       X(1) = -3.0
  458.       X(2) = -1.0
  459.       X(3) =  0.0
  460.       X(4) =  3.0
  461.       X(5) =  4.0
  462.       Y(1) =  7.0
  463.       Y(2) = 11.0
  464.       Y(3) = 26.0
  465.       Y(4) = 56.0
  466.       Y(5) = 29.0
  467.       M = 3
  468.       MM = 2*M
  469.       MM1 = MM - 1
  470.       WRITE(NOUT,15) N,M
  471.    15 FORMAT(5H-N = ,I3,7H    M =,I2)
  472.       CALL QUINAT(N,X,Y,BB,CC,DD,EE,FF)
  473.       DO 20 I = 1,N
  474.          A(I,1) = Y(I)
  475.          A(I,2) = BB(I)
  476.          A(I,3) = CC(I)
  477.          A(I,4) = DD(I)
  478.          A(I,5) = EE(I)
  479.          A(I,6) = FF(I)
  480.    20 CONTINUE
  481.       DO 30 I = 1,MM1
  482.          DIFF(I) = 0.0
  483.          COM(I) = 0.0
  484.    30 CONTINUE
  485.       DO 70 K = 1,N
  486.          DO 35 I = 1,MM
  487.             C(I) = A(K,I)
  488.    35 CONTINUE
  489.       WRITE(NOUT,40) K
  490.    40 FORMAT(40H ---------------------------------------,I3,
  491.      *  45H --------------------------------------------)
  492.       WRITE(NOUT,45) X(K)
  493.    45 FORMAT(F12.8)
  494.       IF (K .EQ. N) WRITE(NOUT,55) C(1)
  495.       IF (K .EQ. N) GO TO 75
  496.       WRITE(NOUT,55) (C(I),I=1,MM)
  497.       DO 50 I = 1,MM1
  498.          IF (ABS(A(K,I)) .GT. COM(I)) COM(I) = ABS(A(K,I))
  499.    50 CONTINUE
  500.    55 FORMAT(6F16.8)
  501.       Z = X(K+1) - X(K)
  502.       DO 60 I = 2,MM
  503.          DO 60 JJ = I,MM
  504.             J = MM + I - JJ
  505.             C(J-1) = C(J)*Z + C(J-1)
  506.    60 CONTINUE
  507.       WRITE(NOUT,55) (C(I),I=1,MM)
  508.       DO 65 I = 1,MM1
  509.          IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 65
  510.          Z = ABS(C(I) - A(K+1,I))
  511.          IF (Z .GT. DIFF(I)) DIFF(I) = Z
  512.    65 CONTINUE
  513.    70 CONTINUE
  514.    75 WRITE(NOUT,80)
  515.    80 FORMAT(41H  MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
  516.       WRITE(NOUT,85) (DIFF(I),I=1,MM1)
  517.    85 FORMAT(5E18.9)
  518.       WRITE(NOUT,90)
  519.    90 FORMAT(42H  MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
  520.       IF (ABS(C(1)) .GT. COM(1)) COM(1) =ABS(C(1))
  521.       WRITE(NOUT,95) (COM(I),I=1,MM1)
  522.    95 FORMAT(5F16.8)
  523.       M = 3
  524.       DO 200 N = 10,100,10
  525.       MM = 2*M
  526.       MM1 = MM - 1
  527.       NM1 = N -1
  528.       DO 100 I = 1,NM1,2
  529.          X(I)   = I
  530.          X(I+1) = I + 1
  531.          Y(I)   = 1.
  532.          Y(I+1) = 0.
  533.   100 CONTINUE
  534.       IF (MOD(N,2) .EQ. 0) GOTO 105
  535.       X(N) = N
  536.       Y(N) = 1.
  537.   105 CONTINUE
  538.       WRITE(NOUT,110) N,M
  539.   110 FORMAT(5H-N = ,I3,7H    M =,I2)
  540.       CALL QUINAT(N,X,Y,BB,CC,DD,EE,FF)
  541.       DO 115 I = 1,N
  542.          A(I,1) = Y(I)
  543.          A(I,2) = BB(I)
  544.          A(I,3) = CC(I)
  545.          A(I,4) = DD(I)
  546.          A(I,5) = EE(I)
  547.          A(I,6) = FF(I)
  548.   115 CONTINUE
  549.       DO 120 I = 1, MM1
  550.          DIFF(I) = 0.0
  551.          COM(I) = 0.0
  552.   120 CONTINUE
  553.       DO 165 K = 1,N
  554.          DO 125 I = 1,MM
  555.             C(I) = A(K,I)
  556.   125    CONTINUE
  557.          IF (N .GT. 10) GOTO 140
  558.          WRITE(NOUT,130) K
  559.   130 FORMAT(40H ---------------------------------------,I3,
  560.      *       45H --------------------------------------------)
  561.          WRITE(NOUT,135) X(K)
  562.   135 FORMAT(F12.8)
  563.          IF (K .EQ. N) WRITE(NOUT,150) C(1)
  564.   140    CONTINUE
  565.          IF (K .EQ. N) GO TO 170
  566.          IF (N .LE. 10) WRITE(NOUT,150) (C(I), I=1,MM)
  567.          DO 145 I = 1,MM1
  568.             IF (ABS(A(K,I)) .GT.  COM(I)) COM(I) = ABS(A(K,I))
  569.   145    CONTINUE
  570.   150 FORMAT(6F16.8)
  571.          Z = X(K+1) - X(K)
  572.          DO 155 I = 2,MM
  573.             DO 155 JJ = I,MM
  574.                J = MM + I - JJ
  575.                C(J-1) = C(J)*Z + C(J-1)
  576.   155    CONTINUE
  577.          IF (N .LE. 10) WRITE(NOUT,150) (C(I), I=1,MM)
  578.          DO 160 I = 1,MM1
  579.             IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 160
  580.             Z = ABS(C(I) - A(K+1,I))
  581.             IF (Z .GT. DIFF(I)) DIFF(I) = Z
  582.   160    CONTINUE
  583.   165 CONTINUE
  584.   170 WRITE(NOUT,175)
  585.   175 FORMAT(41H  MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
  586.       WRITE(NOUT,180) (DIFF(I),I=1,MM1)
  587.   180 FORMAT(5E18.9)
  588.       WRITE(NOUT,185)
  589.   185 FORMAT(42H  MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
  590.       IF (ABS(C(1)) .GT. COM(1)) COM(1) = ABS(C(1))
  591.       WRITE(NOUT,190) (COM(I),I=1,MM1)
  592.   190 FORMAT(5E18.9)
  593.   200 CONTINUE
  594. C
  595. C
  596. C     TEST OF QUINEQ FOLLOWS.
  597. C
  598.       WRITE(NOUT,210)
  599.   210 FORMAT(18H1   TEST OF QUINEQ)
  600.       M = 3
  601.       DO 400 N = 10,100,10
  602.       MM = 2*M
  603.       MM1 = MM - 1
  604.       NM1 = N -1
  605.       DO 300 I = 1,NM1,2
  606.          X(I)   = I
  607.          X(I+1) = I + 1
  608.          Y(I)   = 1.
  609.          Y(I+1) = 0.
  610.   300 CONTINUE
  611.       IF (MOD(N,2) .EQ. 0) GOTO 305
  612.       X(N) = FLOAT(N)
  613.       Y(N) = 1.
  614.   305 CONTINUE
  615.       WRITE(NOUT,310) N,M
  616.   310 FORMAT(5H-N = ,I3,7H    M =,I2)
  617.       CALL QUINEQ(N,Y,BB,CC,DD,EE,FF)
  618.       DO 315 I = 1,N
  619.          A(I,1) = Y(I)
  620.          A(I,2) = BB(I)
  621.          A(I,3) = CC(I)
  622.          A(I,4) = DD(I)
  623.          A(I,5) = EE(I)
  624.          A(I,6) = FF(I)
  625.   315 CONTINUE
  626.       DO 320 I = 1, MM1
  627.          DIFF(I) = 0.0
  628.          COM(I) = 0.0
  629.   320 CONTINUE
  630.       DO 365 K = 1,N
  631.          DO 325 I = 1,MM
  632.             C(I) = A(K,I)
  633.   325    CONTINUE
  634.          IF (N .GT. 10) GOTO 340
  635.          WRITE(NOUT,330) K
  636.   330 FORMAT(40H ---------------------------------------,I3,
  637.      *       45H --------------------------------------------)
  638.          WRITE(NOUT,335) X(K)
  639.   335 FORMAT(F12.8)
  640.          IF (K .EQ. N) WRITE(NOUT,350) C(1)
  641.   340    CONTINUE
  642.          IF (K .EQ. N) GO TO 370
  643.          IF (N .LE. 10) WRITE(NOUT,350) (C(I), I=1,MM)
  644.          DO 345 I = 1,MM1
  645.             IF (ABS(A(K,I)) .GT.  COM(I)) COM(I) = ABS(A(K,I))
  646.   345    CONTINUE
  647.   350 FORMAT(6F16.8)
  648.          Z = 1.
  649.          DO 355 I = 2,MM
  650.             DO 355 JJ = I,MM
  651.                J = MM + I - JJ
  652.                C(J-1) = C(J)*Z + C(J-1)
  653.   355    CONTINUE
  654.          IF (N .LE. 10) WRITE(NOUT,350) (C(I), I=1,MM)
  655.          DO 360 I = 1,MM1
  656.             IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 360
  657.             Z = ABS(C(I) - A(K+1,I))
  658.             IF (Z .GT. DIFF(I)) DIFF(I) = Z
  659.   360    CONTINUE
  660.   365 CONTINUE
  661.   370 WRITE(NOUT,375)
  662.   375 FORMAT(41H  MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
  663.       WRITE(NOUT,380) (DIFF(I),I=1,MM1)
  664.   380 FORMAT(5E18.9)
  665.       WRITE(NOUT,385)
  666.   385 FORMAT(42H  MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
  667.       IF (ABS(C(1)) .GT. COM(1)) COM(1) = ABS(C(1))
  668.       WRITE(NOUT,390) (COM(I),I=1,MM1)
  669.   390 FORMAT(5E18.9)
  670.   400 CONTINUE
  671. C
  672. C
  673. C     TEST OF QUINDF WITH NONEQUIDISTANT KNOTS FOLLOWS.
  674. C
  675.       WRITE(NOUT,410)
  676.   410 FORMAT(50H1         TEST OF QUINDF WITH NONEQUIDISTANT KNOTS)
  677.       N = 5
  678.       X(1) = -3.0
  679.       X(2) = -1.0
  680.       X(3) =  0.0
  681.       X(4) =  3.0
  682.       X(5) =  4.0
  683.       Y(1) =  7.0
  684.       Y(2) = 11.0
  685.       Y(3) = 26.0
  686.       Y(4) = 56.0
  687.       Y(5) = 29.0
  688.       B(1) =  2.0
  689.       B(2) = 15.0
  690.       B(3) = 10.0
  691.       B(4) =-27.0
  692.       B(5) =-30.0
  693.       M = 3
  694.       MM = 2*M
  695.       MM1 = MM - 1
  696.       WRITE(NOUT,415) N,M
  697.   415 FORMAT(5H-N = ,I3,7H    M =,I2)
  698.       CALL QUINDF(N,X,Y,B,CC,DD,EE,FF)
  699.       DO 420 I = 1,N
  700.          A(I,1) = Y(I)
  701.          A(I,2) = B(I)
  702.          A(I,3) = CC(I)
  703.          A(I,4) = DD(I)
  704.          A(I,5) = EE(I)
  705.          A(I,6) = FF(I)
  706.   420 CONTINUE
  707.       DO 430 I = 1,MM1
  708.          DIFF(I) = 0.0
  709.          COM(I) = 0.0
  710.   430 CONTINUE
  711.       DO 470 K = 1,N
  712.          DO 435 I = 1,MM
  713.             C(I) = A(K,I)
  714.   435 CONTINUE
  715.       WRITE(NOUT,440) K
  716.   440 FORMAT(40H ---------------------------------------,I3,
  717.      *  45H --------------------------------------------)
  718.       WRITE(NOUT,445) X(K)
  719.   445 FORMAT(F12.8)
  720.       IF (K .EQ. N) WRITE(NOUT,455) C(1)
  721.       IF (K .EQ. N) GO TO 475
  722.       WRITE(NOUT,455) (C(I),I=1,MM)
  723.       DO 450 I = 1,MM1
  724.          IF (ABS(A(K,I)) .GT. COM(I)) COM(I) = ABS(A(K,I))
  725.   450 CONTINUE
  726.   455 FORMAT(6F16.8)
  727.       Z = X(K+1) - X(K)
  728.       DO 460 I = 2,MM
  729.          DO 460 JJ = I,MM
  730.             J = MM + I - JJ
  731.             C(J-1) = C(J)*Z + C(J-1)
  732.   460 CONTINUE
  733.       WRITE(NOUT,455) (C(I),I=1,MM)
  734.       DO 465 I = 1,MM1
  735.          IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 465
  736.          Z = ABS(C(I) - A(K+1,I))
  737.          IF (Z .GT. DIFF(I)) DIFF(I) = Z
  738.   465 CONTINUE
  739.   470 CONTINUE
  740.   475 WRITE(NOUT,480)
  741.   480 FORMAT(41H  MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
  742.       WRITE(NOUT,485) (DIFF(I),I=1,MM1)
  743.   485 FORMAT(5E18.9)
  744.       WRITE(NOUT,490)
  745.   490 FORMAT(42H  MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
  746.       IF (ABS(C(1)) .GT. COM(1)) COM(1) =ABS(C(1))
  747.       WRITE(NOUT,495) (COM(I),I=1,MM1)
  748.   495 FORMAT(5F16.8)
  749.       M = 3
  750.       DO 600 N = 10,100,10
  751.       MM = 2*M
  752.       MM1 = MM - 1
  753.       P = 0.0
  754.       DO 500 I = 1,N
  755.          P = P + ABS(SIN((FLOAT(I)))) + 0.001*FLOAT(I)
  756.          X(I) = P
  757.          Y(I) = COS((FLOAT(I))) - 0.5
  758.          B(I) = COS((FLOAT(2*I))) - 0.5
  759.   500 CONTINUE
  760.       WRITE(NOUT,510) N,M
  761.   510 FORMAT(5H-N = ,I3,7H    M =,I2)
  762.       CALL QUINDF(N,X,Y,B,CC,DD,EE,FF)
  763.       DO 515 I = 1,N
  764.          A(I,1) = Y(I)
  765.          A(I,2) = B(I)
  766.          A(I,3) = CC(I)
  767.          A(I,4) = DD(I)
  768.          A(I,5) = EE(I)
  769.          A(I,6) = FF(I)
  770.   515 CONTINUE
  771.       DO 520 I = 1, MM1
  772.          DIFF(I) = 0.0
  773.          COM(I) = 0.0
  774.   520 CONTINUE
  775.       DO 565 K = 1,N
  776.          DO 525 I = 1,MM
  777.             C(I) = A(K,I)
  778.   525    CONTINUE
  779.          IF (N .GT. 10) GOTO 540
  780.          WRITE(NOUT,530) K
  781.   530 FORMAT(40H ---------------------------------------,I3,
  782.      *       45H --------------------------------------------)
  783.          WRITE(NOUT,535) X(K)
  784.   535 FORMAT(F12.8)
  785.          IF (K .EQ. N) WRITE(NOUT,550) C(1)
  786.   540    CONTINUE
  787.          IF (K .EQ. N) GO TO 570
  788.          IF (N .LE. 10) WRITE(NOUT,550) (C(I), I=1,MM)
  789.          DO 545 I = 1,MM1
  790.             IF (ABS(A(K,I)) .GT.  COM(I)) COM(I) = ABS(A(K,I))
  791.   545    CONTINUE
  792.   550 FORMAT(6F16.8)
  793.          Z = X(K+1) - X(K)
  794.          DO 555 I = 2,MM
  795.             DO 555 JJ = I,MM
  796.                J = MM + I - JJ
  797.                C(J-1) = C(J)*Z + C(J-1)
  798.   555    CONTINUE
  799.          IF (N .LE. 10) WRITE(NOUT,550) (C(I), I=1,MM)
  800.          DO 560 I = 1,MM1
  801.             IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 560
  802.             Z = ABS(C(I) - A(K+1,I))
  803.             IF (Z .GT. DIFF(I)) DIFF(I) = Z
  804.   560    CONTINUE
  805.   565 CONTINUE
  806.   570 WRITE(NOUT,575)
  807.   575 FORMAT(41H  MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
  808.       WRITE(NOUT,580) (DIFF(I),I=1,MM1)
  809.   580 FORMAT(5E18.9)
  810.       WRITE(NOUT,585)
  811.   585 FORMAT(42H  MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
  812.       IF (ABS(C(1)) .GT. COM(1)) COM(1) = ABS(C(1))
  813.       WRITE(NOUT,590) (COM(I),I=1,MM1)
  814.   590 FORMAT(5E18.9)
  815.   600 CONTINUE
  816. C
  817. C
  818. C     TEST OF QUINAT WITH NONEQUIDISTANT DOUBLE KNOTS FOLLOWS.
  819. C
  820.       WRITE(NOUT,610)
  821.   610 FORMAT(50H1  TEST OF QUINAT WITH NONEQUIDISTANT DOUBLE KNOTS)
  822.       N = 5
  823.       X(1) = -3.
  824.       X(2) = -3.
  825.       X(3) = -1.
  826.       X(4) = -1.
  827.       X(5) = 0.
  828.       X(6) = 0.
  829.       X(7) = 3.
  830.       X(8) = 3.
  831.       X(9) = 4.
  832.       X(10) = 4.
  833.       Y(1) = 7.
  834.       Y(2) = 2.
  835.       Y(3) = 11.
  836.       Y(4) = 15.
  837.       Y(5) = 26.
  838.       Y(6) = 10.
  839.       Y(7) = 56.
  840.       Y(8) = -27.
  841.       Y(9) = 29.
  842.       Y(10) = -30.
  843.       M = 3
  844.       NN = 2*N
  845.       MM = 2*M
  846.       MM1 = MM - 1
  847.       WRITE(NOUT,615) N,M
  848.   615 FORMAT(5H-N = ,I3,7H    M =,I2)
  849.       CALL QUINAT(NN,X,Y,BB,CC,DD,EE,FF)
  850.       DO 620 I = 1,NN
  851.          A(I,1) = Y(I)
  852.          A(I,2) = BB(I)
  853.          A(I,3) = CC(I)
  854.          A(I,4) = DD(I)
  855.          A(I,5) = EE(I)
  856.          A(I,6) = FF(I)
  857.   620 CONTINUE
  858.       DO 630 I = 1,MM1
  859.          DIFF(I) = 0.0
  860.          COM(I) = 0.0
  861.   630 CONTINUE
  862.       DO 670 K = 1,NN
  863.          DO 635 I = 1,MM
  864.             C(I) = A(K,I)
  865.   635 CONTINUE
  866.       WRITE(NOUT,640) K
  867.   640 FORMAT(40H ---------------------------------------,I3,
  868.      *  45H --------------------------------------------)
  869.       WRITE(NOUT,645) X(K)
  870.   645 FORMAT(F12.8)
  871.       IF (K .EQ. NN) WRITE(NOUT,655) C(1)
  872.       IF (K .EQ. NN) GO TO 675
  873.       WRITE(NOUT,655) (C(I),I=1,MM)
  874.       DO 650 I = 1,MM1
  875.          IF (ABS(A(K,I)) .GT. COM(I)) COM(I) = ABS(A(K,I))
  876.   650 CONTINUE
  877.   655 FORMAT(6F16.8)
  878.       Z = X(K+1) - X(K)
  879.       DO 660 I = 2,MM
  880.          DO 660 JJ = I,MM
  881.             J = MM + I - JJ
  882.             C(J-1) = C(J)*Z + C(J-1)
  883.   660 CONTINUE
  884.       WRITE(NOUT,655) (C(I),I=1,MM)
  885.       DO 665 I = 1,MM1
  886.          IF ((K .GE. NN-1) .AND. (I .NE. 1)) GO TO 665
  887.          Z = ABS(C(I) - A(K+1,I))
  888.          IF (Z .GT. DIFF(I)) DIFF(I) = Z
  889.   665 CONTINUE
  890.   670 CONTINUE
  891.   675 WRITE(NOUT,680)
  892.   680 FORMAT(41H  MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
  893.       WRITE(NOUT,685) (DIFF(I),I=1,MM1)
  894.   685 FORMAT(5E18.9)
  895.       WRITE(NOUT,690)
  896.   690 FORMAT(42H  MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
  897.       IF (ABS(C(1)) .GT. COM(1)) COM(1) =ABS(C(1))
  898.       WRITE(NOUT,695) (COM(I),I=1,MM1)
  899.   695 FORMAT(5F16.8)
  900.       M = 3
  901.       DO 800 N = 10,100,10
  902.       NN = 2*N
  903.       MM = 2*M
  904.       MM1 = MM - 1
  905.       P = 0.0
  906.       DO 700 I = 1,N
  907.          P = P + ABS(SIN(FLOAT(I)))
  908.          X(2*I-1) = P
  909.          X(2*I)   = P
  910.          Y(2*I-1) = COS(FLOAT(I)) - 0.5
  911.          Y(2*I)   = COS(FLOAT(2*I)) - 0.5
  912.   700 CONTINUE
  913.       WRITE(NOUT,710) N,M
  914.   710 FORMAT(5H-N = ,I3,7H    M =,I2)
  915.       CALL QUINAT(NN,X,Y,BB,CC,DD,EE,FF)
  916.       DO 715 I = 1,NN
  917.          A(I,1) = Y(I)
  918.          A(I,2) = BB(I)
  919.          A(I,3) = CC(I)
  920.          A(I,4) = DD(I)
  921.          A(I,5) = EE(I)
  922.          A(I,6) = FF(I)
  923.   715 CONTINUE
  924.       DO 720 I = 1, MM1
  925.          DIFF(I) = 0.0
  926.          COM(I) = 0.0
  927.   720 CONTINUE
  928.       DO 765 K = 1,NN
  929.          DO 725 I = 1,MM
  930.             C(I) = A(K,I)
  931.   725    CONTINUE
  932.          IF (N .GT. 10)  GOTO 740
  933.          WRITE(NOUT,730) K
  934.   730 FORMAT(40H ---------------------------------------,I3,
  935.      *       45H --------------------------------------------)
  936.          WRITE(NOUT,735) X(K)
  937.   735 FORMAT(F12.8)
  938.          IF (K .EQ. NN) WRITE(NOUT,750) C(1)
  939.   740    CONTINUE
  940.          IF (K .EQ. NN) GO TO 770
  941.          IF (N .LE. 10) WRITE(NOUT,750) (C(I), I=1,MM)
  942.          DO 745 I = 1,MM1
  943.             IF (ABS(A(K,I)) .GT.  COM(I)) COM(I) = ABS(A(K,I))
  944.   745    CONTINUE
  945.   750 FORMAT(6F16.8)
  946.          Z = X(K+1) - X(K)
  947.          DO 755 I = 2,MM
  948.             DO 755 JJ = I,MM
  949.                J = MM + I - JJ
  950.                C(J-1) = C(J)*Z + C(J-1)
  951.   755    CONTINUE
  952.          IF (N .LE. 10) WRITE(NOUT,750) (C(I), I=1,MM)
  953.          DO 760 I = 1,MM1
  954.             IF ((K .GE. NN-1) .AND. (I .NE. 1)) GO TO 760
  955.             Z = ABS(C(I) - A(K+1,I))
  956.             IF (Z .GT. DIFF(I)) DIFF(I) = Z
  957.   760    CONTINUE
  958.   765 CONTINUE
  959.   770 WRITE(NOUT,775)
  960.   775 FORMAT(41H  MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
  961.       WRITE(NOUT,780) (DIFF(I),I=1,MM1)
  962.   780 FORMAT(5E18.9)
  963.       WRITE(NOUT,785)
  964.   785 FORMAT(42H  MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
  965.       IF (ABS(C(1)) .GT. COM(1)) COM(1) = ABS(C(1))
  966.       WRITE(NOUT,790) (COM(I),I=1,MM1)
  967.   790 FORMAT(5E18.9)
  968.   800 CONTINUE
  969. C
  970. C
  971. C     TEST OF QUINAT WITH NONEQUIDISTANT KNOTS, ONE DOUBLE KNOT,
  972. C        ONE TRIPLE KNOT, FOLLOWS.
  973. C
  974.       WRITE(NOUT,805)
  975.       WRITE(NOUT,810)
  976.   805 FORMAT(51H1         TEST OF QUINAT WITH NONEQUIDISTANT KNOTS,)
  977.   810 FORMAT(40H             ONE DOUBLE, ONE TRIPLE KNOT)
  978.       N = 8
  979.       X(1) = -3.
  980.       X(2) = -1.
  981.       X(3) =  -1.
  982.       X(4) =  0.
  983.       X(5) =  3.
  984.       X(6) = 3.
  985.       X(7) = 3.
  986.       X(8) = 4.
  987.       Y(1) =  7.
  988.       Y(2) = 11.
  989.       Y(3) = 15.
  990.       Y(4) = 26.
  991.       Y(5) = 56.
  992.       Y(6) = -30.
  993.       Y(7) =  -7.
  994.       Y(8) =  29.
  995.       M = 3
  996.       MM = 2*M
  997.       MM1 = MM - 1
  998.       WRITE(NOUT,815) N,M
  999.   815 FORMAT(5H-N = ,I3,7H    M =,I2)
  1000.       CALL QUINAT(N,X,Y,BB,CC,DD,EE,FF)
  1001.       DO 820 I = 1,N
  1002.          A(I,1) = Y(I)
  1003.          A(I,2) = BB(I)
  1004.          A(I,3) = CC(I)
  1005.          A(I,4) = DD(I)
  1006.          A(I,5) = EE(I)
  1007.          A(I,6) = FF(I)
  1008.   820 CONTINUE
  1009.       DO 830 I = 1,MM1
  1010.          DIFF(I) = 0.0
  1011.          COM(I) = 0.0
  1012.   830 CONTINUE
  1013.       DO 870 K = 1,N
  1014.          DO 835 I = 1,MM
  1015.             C(I) = A(K,I)
  1016.   835 CONTINUE
  1017.       WRITE(NOUT,840) K
  1018.   840 FORMAT(40H ---------------------------------------,I3,
  1019.      *  45H --------------------------------------------)
  1020.       WRITE(NOUT,845) X(K)
  1021.   845 FORMAT(F12.8)
  1022.       IF (K .EQ. N) WRITE(NOUT,855) C(1)
  1023.       IF (K .EQ. N) GO TO 875
  1024.       WRITE(NOUT,855) (C(I),I=1,MM)
  1025.       DO 850 I = 1,MM1
  1026.          IF (ABS(A(K,I)) .GT. COM(I)) COM(I) = ABS(A(K,I))
  1027.   850 CONTINUE
  1028.   855 FORMAT(6F16.8)
  1029.       Z = X(K+1) - X(K)
  1030.       DO 860 I = 2,MM
  1031.          DO 860 JJ = I,MM
  1032.             J = MM + I - JJ
  1033.             C(J-1) = C(J)*Z + C(J-1)
  1034.   860 CONTINUE
  1035.       WRITE(NOUT,855) (C(I),I=1,MM)
  1036.       DO 865 I = 1,MM1
  1037.          IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 865
  1038.          Z = ABS(C(I) - A(K+1,I))
  1039.          IF (Z .GT. DIFF(I)) DIFF(I) = Z
  1040.   865 CONTINUE
  1041.   870 CONTINUE
  1042.   875 WRITE(NOUT,880)
  1043.   880 FORMAT(41H  MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
  1044.       WRITE(NOUT,885) (DIFF(I),I=1,MM1)
  1045.   885 FORMAT(5E18.9)
  1046.       WRITE(NOUT,890)
  1047.   890 FORMAT(42H  MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
  1048.       IF (ABS(C(1)) .GT. COM(1)) COM(1) =ABS(C(1))
  1049.       WRITE(NOUT,895) (COM(I),I=1,MM1)
  1050.   895 FORMAT(5F16.8)
  1051. C
  1052. C
  1053. C     TEST OF QUINAT WITH NONEQUIDISTANT KNOTS, TWO DOUBLE KNOTS,
  1054. C        ONE TRIPLE KNOT,FOLLOWS.
  1055. C
  1056.       WRITE(NOUT,905)
  1057.       WRITE(NOUT,910)
  1058.   905 FORMAT(51H1         TEST OF QUINAT WITH NONEQUIDISTANT KNOTS,)
  1059.   910 FORMAT(40H             TWO DOUBLE, ONE TRIPLE KNOT)
  1060.       N = 10
  1061.       X(1) = 0.
  1062.       X(2) = 2.
  1063.       X(3) = 2.
  1064.       X(4) = 3.
  1065.       X(5) = 3.
  1066.       X(6) = 3.
  1067.       X(7) = 5.
  1068.       X(8) = 8.
  1069.       X(9) = 9.
  1070.       X(10)= 9.
  1071.       Y(1) = 163.
  1072.       Y(2) = 237.
  1073.       Y(3) = -127.
  1074.       Y(4) = 119.
  1075.       Y(5) = -65.
  1076.       Y(6) = 192.
  1077.       Y(7) = 293.
  1078.       Y(8) =  326.
  1079.       Y(9) = 0.
  1080.       Y(10)= -414.0
  1081.       M = 3
  1082.       MM = 2*M
  1083.       MM1 = MM - 1
  1084.       WRITE(NOUT,915) N,M
  1085.   915 FORMAT(5H-N = ,I3,7H    M =,I2)
  1086.       CALL QUINAT(N,X,Y,BB,CC,DD,EE,FF)
  1087.       DO 920 I = 1,N
  1088.          A(I,1) = Y(I)
  1089.          A(I,2) = BB(I)
  1090.          A(I,3) = CC(I)
  1091.          A(I,4) = DD(I)
  1092.          A(I,5) = EE(I)
  1093.          A(I,6) = FF(I)
  1094.   920 CONTINUE
  1095.       DO 930 I = 1,MM1
  1096.          DIFF(I) = 0.0
  1097.          COM(I) = 0.0
  1098.   930 CONTINUE
  1099.       DO 970 K = 1,N
  1100.          DO 935 I = 1,MM
  1101.             C(I) = A(K,I)
  1102.   935 CONTINUE
  1103.       WRITE(NOUT,940) K
  1104.   940 FORMAT(40H ---------------------------------------,I3,
  1105.      *  45H --------------------------------------------)
  1106.       WRITE(NOUT,945) X(K)
  1107.   945 FORMAT(F12.8)
  1108.       IF (K .EQ. N) WRITE(NOUT,955) C(1)
  1109.       IF (K .EQ. N) GO TO 975
  1110.       WRITE(NOUT,955) (C(I),I=1,MM)
  1111.       DO 950 I = 1,MM1
  1112.          IF (ABS(A(K,I)) .GT. COM(I)) COM(I) = ABS(A(K,I))
  1113.   950 CONTINUE
  1114.   955 FORMAT(6F16.8)
  1115.       Z = X(K+1) - X(K)
  1116.       DO 960 I = 2,MM
  1117.          DO 960 JJ = I,MM
  1118.             J = MM + I - JJ
  1119.             C(J-1) = C(J)*Z + C(J-1)
  1120.   960 CONTINUE
  1121.       WRITE(NOUT,955) (C(I),I=1,MM)
  1122.       DO 965 I = 1,MM1
  1123.          IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 965
  1124.          Z = ABS(C(I) - A(K+1,I))
  1125.          IF (Z .GT. DIFF(I)) DIFF(I) = Z
  1126.   965 CONTINUE
  1127.   970 CONTINUE
  1128.   975 WRITE(NOUT,980)
  1129.   980 FORMAT(41H  MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
  1130.       WRITE(NOUT,985) (DIFF(I),I=1,MM1)
  1131.   985 FORMAT(5E18.9)
  1132.       WRITE(NOUT,990)
  1133.   990 FORMAT(42H  MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
  1134.       IF (ABS(C(1)) .GT. COM(1)) COM(1) =ABS(C(1))
  1135.       WRITE(NOUT,995) (COM(I),I=1,MM1)
  1136.   995 FORMAT(5F16.8)
  1137.       STOP
  1138.       END
  1139.